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Abstract 

»vj ■ We apply two different monotonically convergent optimization algo- 

rithms to the control of molecular rotational dynamics by laser pulses. 
This example represents a quantum control problem where the interaction 

1-^ ■ of the system with the external field is non-linear. We test the validity 

and accuracy of the two methods on the key control targets of producing 
molecular orientation and planar delocalization at zero temperature, and 
maximizing permanent alignment at non-zero temperature. 



_1j' 1 Introduction 



Optimal control theory is nowadays a mature mathematical discipline with a 
wide range of applications in science and engineering [XJ- The technique has 
^^ . been used with success in quantum mechanics since the beginning of the 1990s 

l/~j ' [21 [31 131 [S] to control spins, atoms and molecules by external electromagnetic 

^O ' fields. Control problems can be tackled by two different types of approaches, 

t:;;^.' ; geometric [^ [3 [H] and numerical methods [3 [TOl [H [H [H [HI [IS] for quantum 

f^ . systems of low and high dimension, respectively. It is this second aspect which 

Cn ' is at the core of this article. Numerical optimal control algorithms can roughly 

be divided into Gradient ascent algorithms [Hj and Krotov's method [TUl [HI 
[2U1 [2T1 [T71 [TSl FT^ . The latter guarantees monotonic convergence independent 
of the specific choice of optimization functional, type of interaction between 
system and external control, and equations of motion. In the quantum control 
C^ , literature, Krotov's method was first established for dipole transitions, where 

the interaction of the system with the control field is linear [101 HIl HH US] ■ In 
recent years, several modifications to the known algorithms have been brought 
forward to account for the non-linear case, a problem which arises naturally in 
a variety of control problems in atomic and molecular physics. In particular it 
occurs when the iirtensity of the laser field is sufficieirtly large, so that the linear 
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model is no longer a good approximation of the dynamical system. While the 
generalization is straightforward for gradient algorithms, the extension of the 
monotonia approach is more involved |24[ [23 [TH] . Here, our goal is to explore 
the efficiency of two different schemes of monotonically convergent optimization 
algorithms for the control of a molecule interacting non-linearly with the control 
field. 

The control of molecular rotation [SHIlTlllHllMllSnilSIllMllSSl, for which 
such non- linear models are well established |34[|35) . is used as a testbed case to 
analyze the features of these algorithms. A first modification of a monotonically 
convergent algorithm to account for a non-linear interaction with the control 
assumes the cost to be quadratic in the field and decomposes the control into n 
components for a nonlinearity of order n |24| . The decomposition leads to 2n 
Schrodinger equations that need to be solved, n for the wave function and n for 
the adjoint state. This can be numerically costly. The approach was successfully 
applied to the control of molecular orientation and alignment [36J . At the same 
time, some of us proposed a new algorithm using only one component of the wave 
function [33]. This comes at the price of changing the cost functional. Instead of 
penalizing the intensity of the field, i.e., the square of the control parameter, it 
penalizes a higher exponent, the value of which depends on the order of the non- 
linearity. The algorithms of Refs. [211 [2S] have recently been compared [37]. In 
the case of a two-color control strategy for molecular orientation, it was shown 
that the efficiency of the two optimized solutions designed by the two algorithms 
was similar. In parallel, it has been mentioned that the Krotov method allows 
for constructing a monotonically convergent algorithm with the standard cost 
functional penalizing the field intensity |19| . Here, we examine this claim and 
perform an extensive comparison with the Lapert algorithm |25| . We analyze 
the efficiency, numerical cost and structure of the optimized solutions obtained 
by the two approaches. The rotational dynamics of a diatomic molecule driven 
by an electromagnetic field will be used as an illustrative example. 

The remainder of this paper is organized as follows. The molecular model 
is presented in Sec. [21 Section [31 is devoted to the application of the two 
algorithms to the control objectives of controlling molecular orientation and 
planar delocalization at zero temperature, and producing permanent alignment 
at non-zero temperature. We conclude in Sec.[H Appendix \X\ summarizes briefly 
the two optimization algorithms. 

2 The model 

We consider the control of the rotational dynamics of the linear CO molecule 
described in a rigid rotor approximation and driven by the electric field E{t). 
The field is expressed as follows: 

E(t) = exit)coa{LL>t + ^x)ex + ez{t)cos{ujt + ^y)ey + 

e^(t)cos(a;i + $z)e;, (1) 

where ei,(t), e^, and <I>,y, v = x, y^ z are the amplitude, the unit vector and 
the phase along the i^-axis, respectively. The Hamiltonian of the system can be 
written as ^., ,35J: 

H = BJ'^ + /i • E{t) + a-E^{t)+f3- E^{t), (2) 



B (cm-^) 


Mo (a.u.) 


a|| (a.u.) 


a± (a.u.) 


/3|l (a.u.) 


(3± (a.u.) 


1.9312 


0.112 


15.65 


11.73 


28.35 


6.64 



Table 1: Numerical values of the different molecular parameters. 



where B is the rotational constant. The first term of the right-hand side of Eq. 
is the field-free rigid-rotor Hamiltonian. Its eigenstates are the spherical 
harmonics denoted by \j,m), with j > and \m\ < j. The operators /i, a and 
/3 are associated, respectively, to the permanent dipole moment and the polar- 
izability and hyperpolarizability tensors. The spatial position of the diatomic 
molecule is given in the laboratory frame by the spherical coordinates (9, (p). 

We first study the interaction of the molecule with a non-resonant laser 
field, polarized linearly along the z-axis of the laboratory frame. In this case, 
the variable 9 is the angle between the molecular axis and the polarization vector 
of the electric field. The Hamiltonian 1^ then simplifies to 



H 



BJ'^ - fio cos 9E^ (t) - - [ Aa cos^ 9 + a±l] El (t) 



3Pa_) cos^ 9 + 3/3 !_ cos 9] El{t), 



(3) 



where Aa = aii — a^. The parameter /io is the permanent dipole moment and 
the coefficients an. aj_, Pn and I3±_ denote, respectively, the polarizability and 
hyperpolarizability components of the molecule with the labels || and _L indi- 
cating the components parallel and perpendicular to the internuclear axis. The 
numerical values used in our simulations for the different molecular parameters 
are reported in Tab. [TJ For details see Ref. [25]. If we further assume that the 
frequency uj of the laser field is much higher than the rotational frequencies and 
non-resonant with respect to all rovibronic transitions, we can average over the 
fast oscillations of the electric field in Eq. ([3)) and obtain |35l [M] : 



H{t) = BJ"^ -^[/^acos^ 9 + ai_t\el{t) 



-- [(/3|| - 3P ^) cos^ 9 + 3P^ cos 9] el{t). 



(4) 



As a second example we consider the interaction of the CO molecule with 
a pulse that is elliptically polarized in the [x, j/)-plane. We neglect here the 
hyperpolarizability term of the interaction since it does not play a quantitative 
role in this case. After optical-cycle averaging as above, the corresponding 
Hamiltonian is expressed as: 



H 



BJ^ - - [{l\a cos^ 9, J, + a ^t) €l{t)+ 
(Aacc 



? 9y + a^^) €l{t)+ 
2Aacos($j; — ^y) cos6xCos9y€x{t)ey{t)] , 



(5) 



where cos 9x = sin 6 cos Lp and cos 9y — sin 6 sin Lp. 

In the case of zero rotational temperature (T = K) . the time evolution of 
the system is described by the time-dependent Schrodinger equation, 



z||V(t))-i/(i)|^(i)), 



(6) 



where \ip{t)) is the wave function of the system at time t. The Liouville equation 
is used to describe the time evolution for T ^ K: 

^§-^pit) = [Hit),pit)], (7) 

where p{t) represents the density matrix associated with the system at time t. 
Equations © and ([7]) are solved numerically with the split operator algorithm 
|38| . The Hamiltonian is represented in spherical harmonics \j,m) where all 
matrix elements are known analytically. We use atomic units unless otherwise 
specified. 



3 Numerical results 

We explore three different control targets presenting a comparative study of 
the Krotov and Lapert algorithms. The technical details of the algorithms are 
briefly reviewed in Appendix |^ 

3.1 Orientation dynamics driven by a linearly polarized 
field 

We first investigate the control of molecular orientation by a field linearly po- 
larized along the z-axis. In this case, the dynamics is described by Eq. ([3]). The 
control duration tf is chosen to be equal to one rotational period Tp^r of the 
molecule, Tpcr ~ 8.6 ps. We consider a finite Hilbert space of size jmax — 15, 
which is sufficient for the intensity of the laser field used here. The expectation 
value (cos 9) is taken as a quantitative measure of the orientation. The molecule 
is oriented when | (cos 0)1 ~ 1. Following Refs. [31130], we do not maximize this 
expectation value but a target state 1-0/) maximizing |(cos0)| in a sub- Hilbert 
space of finite dimension defined by j < jf. The details of \tpf) can be found 
in Refs. |39[ I40j. Figure [T] shows the projection of the target state onto the 
eigenstates \j, 0) of the molecule, with jf = 4. A Gaussian pulse of 144 fs full 
width at half maximum (FWHM), centered at to = Tper/5 is taken as guess 
field for all optimizations discussed in this section: 

E{t)^E^e-^^, (8) 

where ^eQcE^ — 10^^ W/cm^ is the peak intensity of the laser field and the 
parameter a is defined such that FWHM = 2v2hi2(T. This choice of guess field 
is standard in the control of molecular orientation |39| . It provides an efficient 
initial solution with a population transfer from the state |0, 0) to a superposition 
of |j, 0) states with j = 0, 1, • • • , j/. The cost functional is defined by: 

C = |(^/|^(i/))P - A / \E{t) ~ Eref{t)r/Sit)dt, (9) 

Jo 

where n — 2 and n = 4 for the Krotov and Lapert algorithms, respectively. The 
parameter A penalizes the pulse energy and Eref{t) is a pulse reference. The 
function 5*, given by S{t) — sin (Trt/tf), suppresses pulse amplitude at the be- 
ginning and end of the time window, ensuring a pulse that is smoothly switched 
on and off. We denote the final fidelity of the control by Ct^ = |(V'/IV'(^/))P- 




Figure 1: (color online) Population of the target state maximizing orientation 
of the CO molecule along the z-axis. The corresponding wave function is given 
by IV'/) « 0.34|0, 0) + 0.54|1, 0) + 0.56|2, 0) + 0.46|3, 0) + 0.25|4, 0) 



We first analyze the role of the parameter A in the two algorithms. The 
results reported in Tab. [5] show that the Krotov algorithm requires smaller A 
values to converge to a high fidelity Ct^ than the Lapert method. In order to 
observe a convergence with realistic optimized pulses, the parameter A should 
be larger than 10^ for the Lapert method and lower than 0.1 for the Krotov 
one. This difference is easily understood from the fact that, in the Krotov 
formulation, the running cost is a quadratic function of the electric field, cf. 
Eq. ([9]), while in the algorithm by Lapert et al., this power is 4. Since in atomic 
units, \E{t) — E-cci\ < 1, the same order of magnitude for the two running costs 
is obtained for A'^p > lO^A^'. Note that the parameter A^'' has to satisfy 
the relation (23) in order for the Krotov algorithm to be monotonic. The fact 
that there is no constraint in the choice of A'^p allows more flexibility in the 
use of the Lapert algorithm. However, one should keep in mind that large 
values of A'^p are required in order to avoid fast oscillations in the optimized 
solution. In this section, the two coefficients will be fixed at A'^p = 5 x 10^ and 
A^' = 5 X 10^^. The evolution of the final fidelity Ctj, as a function of the 
number of iterations and the CPU time is displayed in Fig. [21 While the Lapert 
formulation converges faster initially, the insets of Fig. [5] show that the Krotov 
algorithm becomes faster when the fidelity is close to 1. 



A 




cr, 


pip 


pkr 

-^max 


510'^ 


0.9892 


0.0205 


5.510-3 


5.510-3 


sio-* 


0.9993 


0.0205 


5.510-3 


10-2 


5102 


0.9996 


0.0205 


5.510~3 


2.810-2 


510"! 


— 


0.5789 


— 


5.510-3 


510-2 


— 


0.9959 


— 


5.510-3 


510-3 


— 


0.9944 


— 


5.510-3 



Table 2: Analysis of the convergence of the two algorithms with respect to the 
parameter A. C^ and Cj" are the fidelities obtained from Krotov and Lapert 
algorithms, respectively. E]^^^ and -E'^'^x correspond to the maximum amplitude 
in absolute value of the two optimized solutions. The number of iterations is set 
to 20. In the Lapert algorithm, small values of A lead to very fast oscillations of 
the corresponding optimized field, which are not physically relevant. Therefore 
results are not indicated in columns 2 and 4 when A is smaller than 5 x 102. 



The Lapert algorithm is however more costly in terms of computer time 
(CPU time). The faster convergence of Krotov's method in terms of CPU 
time is not surprising since in the Lapert formulation, at each iteration, the 
roots of a polynomial of power 3 need to be determined in order to compute 
the updated new pulse, see the appendix IA.2I for details. Figure [3] compares 
the optimized pulses obtained from the Lapert algorithm. Fig. [3^ and from 
the Krotov one. Fig. [3J;. Figure [SJd displays the guess field considered for the 
two algorithms. Note that the structure of the Krotov solution is very simple, 
since the optimized field is mainly composed of the guess pulse plus additional 
small deformations. The solution designed by the Lapert algorithm is rather 
more complex, in the sense that fast oscillations appear between the middle and 
the end of the optimization time interval. As already pointed out in [25J , this 
behavior seems to be quite general with a cost functional penalizing the power 
4 of the field. Note that spectral filters can be added to avoid such oscillatory 
structures fJTlH^ . 

The dynamics induced by the two optimized pulses is plotted in Fig. U) 
More precisely, figure |3] displays the projection of the wave function onto the 
molecular eigenstates as a function of time. The two dynamics show similar 
features in the first fifth of the optimization time interval, [0,0.2 x tf]. Most 
of the population remains in the ground state, \j = 0,m = 0), since, in this 
time interval, both optimized pulses are almost zero. At time t = 0.2 x tf, 
both optimized solutions contain a kick which leads to a superposition of states 
with j = 0, . . . , 3. For the dynamics induced by the Krotov optimized pulse, 
we observe that most of the population is concentrated in states |j = 3, 0) and 
\j = 4,0) during the time interval [0.4 x t/,0.6 x tf], while these states are 
populated for t £ [0.6 x t/,0.8 x tf] in the Lapert case. In the time interval 
[0.6 X tf, 0.8 X tf], the dynamics induced by the Krotov optimized field shows a 
superposition of states |1,0), |2,0) and |3,0), more or less similar to the target 
state. In the final step, i.e., in the time interval [O.Sxi/, tf], the small oscillations 
of the Lapert and Krotov fields are responsible for the complete transfer of the 
superpositions to the target state. 



u 



-^ 0.6 








12 15 18 21 24 27 30 
iterations 




1 



2 3 4 5 6 7 
CPU time (min) 



10 



Figure 2: (color online) Convergence of the Lapert and Krotov algorithms mea- 
sured by the final cost Ct^ plotted as a function of the number of iterations 
(top panel) and Ctf plotted as a function of the CPU time (bottom panel). 
The results for the Lapert formulation are depicted in black solid line with cir- 
cles, while a blue dashed line with stars represents the results with the Krotov 
algorithm. 



3.2 Delocalization in the {x,y)- plane 

The second example is dedicated to controlling the orientation of the angular 
momentum of the CO molecule along the z-axis of the laboratory frame. The 
degree of orientation is given here by (Jz)/ \/{J^), where J^ is the z- compo- 
nent of the angular momentum J. This aspect has been recently studied in a 
series of works, both theoretically [151 HH H5] and experimentally [^ ITfl US] . 
In particular, it has been shown in Ref. |45j that orientation of the angular 
momentum can be achieved by a sequence of two short laser pulses, properly 
delayed and polarized at 45 degrees with respect to each other. Here we revisit 
this control problem using the two monotonically convergent algorithms. An el- 
liptical polarization is considered to realize this orientation. The corresponding 
Hamiltonian is given by Eq. ^. Since the angular momentum of a diatomic 
molecule is classically orthogonal to the molecular axis, its alignment along the 
laboratory z- axis is equivalent to a delocalization of the molecular axis in the 
(x, y)- plane. This delocalization can also be interpreted as a minimization of 
the expectation value (cos^ 9) [H] . 

Let us consider a state of the form |j, ±j). Straightforward computation 
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Figure 3: (color online) Comparison of the optimized pulses obtained from the 
Krotov (c) and Lapert (a) algorithms. The panel (b) depicts the initial field 
used for the two algorithms. 



shows that 



(j,±j|cos2 0|j,±j)- 



2j 



(10) 



When j ^ oo, the right hand side of Eq. pUj) converges to its minimum value, 0. 
Therefore, in a sub-Hilbert space of finite dimension, the states |j, ±j) minimize 
(cos^ 6) for large j. In other words, these states maximize the delocalization of 
the molecular axis in the (a;,y)-plane. Consequently, the states \j,j) and the 
states \j,—j) maximize and minimize (Jz), respectively f45J. 

At T = K, the initial state is |0,0) and |4,4) is taken as the target state. 
The minimum expectation value that can be reached with this choice is of the 
order of (4, ±4| cos^ 6*14, ±4) « 0.1. Note that the more we increase j, the 
better the delocalization becomes. However, both difhculty and numerical cost 
increase with j. Therefore, target states of the form |j, ±j) with j > 4 will not 
be analyzed. 

When the relative phase ^^ — ^y is chosen so that the cross term of the 
Hamiltonian (O vanishes, the dynamics cannot distinguish the states \j,j) and 
\jj~j) (see [44J for the analytical proof). Here, in order to get a completely 
controllable system, the relative phase ^^ — ^y is set to 7r/4. The guess field is 
constructed as a series of Gaussian pulses of 150 fs FWHM for each component, 
ex and ey. We have chosen X^^ — 10^ and A^"^ =0.1. Our choice of a large A 
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Figure 4: (color online) Projection of the time-dependent wave function onto 
the eigenstates of the system. The top and bottom panels represent respectively 
the dynamics induced by the pulses optimized with the Krotov and Lapert 
algorithms. 

value for the Lapert algorithm is motivated by the fact that small values induce 
very fast oscillations of the optimized control field, which are physically and 
numerically not very interesting. For 50 iterations, the target state is reached 
with a probability of the order of 0.99 for the two algorithms. The convergence 
of the cost Ctf (top panel), plotted as a function of the number of iterations, is 
shown in Fig. [5] for the two algorithms. The corresponding optimized pulses are 
shown in the bottom panels. The convergence to the target state is similar to the 
one observed in the first example. In particular, the Lapert algorithm converges 
faster initially, the Krotov algorithm remaining more efficient, specifically in 
terms of CPU time. In addition. Fig. [5] clearly shows that the structure of the 
Krotov solution is simpler than the Lapert one. 



3.3 Control of a thermal rotational sample 

The last example concerns the control of a thermal sample at non-zero temper- 
ature. The initial state is given by the Boltzmann distribution, which can be 
written as follows: 



Po 



1 ^ 



] m= 



'\j,m){j,m\ 



(11) 



(a) 
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Figure 5: (Color online) (a): Convergence of the Lapert and the Krotov al- 
gorithms measured by the final cost Ctj,- Ctj — |(4,4|'(/;(i/))p (top panel) is 
plotted as a function of the number of iterations, (b): x- components of the 
optimized fields, (c): y- components the optimized fields. The black and blue 
(dark gray) curves depict respectively the results derived with the Lapert and 
the Krotov algorithms. 



where T is the temperature, ks the Boltzmann constant and Z the partition 
function, which is expressed as 



^-E 



3 



Bj(j + 1) 



(12) 



The control aims at reaching the state that maximizes the permanent alignment 
of the molecule along the z-axis. The molecular alignment is measured by the 
expectation value (cos^ 6), which can be written as the sum of two terms: 



(cos 



(cos 



(cos 



(13) 

m.j' ra'^^ jm,j' m' 



where (COS^ 0)^ = Ej,m,m' Pj^Jm'CjmJm' and (C0S2 0)c = 'Ej^j',m,m' Pjr 

The coefficients Cjmj'm' denote the matrix elements of the operator cos^ 6. Par- 
titioning the alignment measured into diagonal and off-diagonal terms (with re- 
spect to the quantum number j) reveals interesting physical information about 
the rotational dynamics. While (cos^ 6)p provides a direct measure of the rota- 
tional population, (cos^ 9)c leads to the temporal evolution of the coherences. 
By definition, (cos^ 9)p is constant when the pulse is switched off. In some 
applications, it can be interesting to maximize only the permanent alignment. 
For this purpose, we use the strategy proposed in Ref. ^\. Considering a 
sub-Hilbert space Hj, of finite dimension defined by the condition j < jf, we 
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introduce the diagonal projection of the cos^ 9 operator: 

cos^e.p= Y^ \j,m){j,m\cos^O\j,m'){j,m'\, (14) 

such that (cos^ Op) = (cos^ 9)p. The target state pj of the control problem 
is defined as the density matrix maximizing (cos^ 9p) and reachable from the 
initial state po- Due to the constraint of unitary evolution, the density matrices 
Po and pf have the same spectrum. In addition, it can be shown that, in this 
subspace, the two operators cos^ Op and pf can be simultaneously diagonalized. 
One therefore deduces that 

N 

max[(cos^ Op)] = ^ Xk^k, (15) 

fc=i 

where xi ^ X2 1^ • • ■ < Xn and ui < lj2 < ■ ■ ■ < lon are the eigenvalues of 
cos^ Op and pf , respectively. The integer N is the dimension of TLj^ . If we denote 
by \xk) the eigenvectors of cos^ 9p, pf becomes 

N 
k=l 

Since the subspaces of a given parity of j are not coupled by the operators 
cos^6'2.„r, the subdivision Hi, = '}{Y""'^"> q) 'Hy has to be considered to 
properly define the target state, see Ref. ^^ for details of this construction. 

Figure [5] displays the partial trace Trj[-] with respect to j of the target state 
Pf, with jf = 4, for the CO molecule at T = 5K. For a given value m, this 
trace is defined by X^i"?™! P'jmjm- ^or comparison, we have also plotted the 
same distribution for the initial state po- One clearly sees in Fig. [6] that the 
optimal distribution is narrowed compared to the thermal one. Since pf is a 
diagonal matrix, there is no coherence and we get (cos^)c = 0. At T = 
K, the maximum permanent alignment is equal to 0.6. This maximum is a 
temperature-dependent function, and for T 7^ K, max[(cos^ Op)] < 0.6. For 
example at T = 5 K, max[(cos^ Op)] = 0.518. 

We use Eq. (O and the same guess pulse as in Sec. 13.21 The parameter jf is 
fixed to 4. We have chosen A = 2 x 10~^ and 5 x 10^ for the Krotov and Lapert 
algorithms, respectively. Figure [7] compares the permanent alignment dynamics 
for the two algorithms at T = 5 K. The pulse is switched off at tf = Tp^r- 
The dynamics are found to be step-like, such that (cos^ 0)p is either constant or 
varies suddenly. A permanent alignment of the order of 0.47 and 0.49 is reached 
for the Lapert and the Krotov algorithms, respectively for 500 iterations. In this 
example, the comparison of the convergence measured by the final cost is not 
shown. We observe as in Fig. [5] the same behavior for the Lapert and Krotov 
algorithms. As illustrated in Fig. [7J the optimized solution designed by the 
Krotov algorithm is simpler than the one given by the Lapert algorithm. 

4 Summary 

For key control problems of the rotational dynamics, we have designed in this 
work different optimized pulses from two monotonically convergent algorithms. 
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-4-3-2-1 1 2 3 4 

m 

Figure 6: (color online) Partial trace Trj[-] with respect to the quantum number 
j of initial and target state distributions at T = 5 K. The target distribution 
corresponds to pf , (red or dark gray bar) which is the density matrix maximizing 
the permanent alignment at T = 5 K. The partial trace of the initial density 
matrix is in gray. The parameter jf is set to 4. 

the Lapert and the Krotov ones. Our numerical findings confirmed by the three 
examples discussed in this work are as follows: 

1. The final fidelities reached by the two algorithms are very similar. 

2. The Lapert algorithm is somehow more fiexible compared to the Kro- 
tov one in the sense that the parameter A can be chosen without any 
constraints, while this parameter has to satisfy the relation (23) hi the 
Krotov approach in order to ensure monotonic convergence. 

3. The Krotov algorithm is more efficient than the Lapert one in terms of 
CPU time. 

4. The optimized Krotov field has a simpler structure than the Lapert one, 
which generally presents some unwanted oscillatory behaviors. 

Having in mind the work of Ref. [49J , an open question is now the application 
of these optimal control algorithms to more complex systems. The computation 
of the optimal field allowing the cooling of rovibrational dynamics could be an 
interesting test case, in particular because non unitary processes have to be 
taken into account. 
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Figure 7: (Color online) (a): Permanent alignment dynamics during and af- 
ter the application of the optimized pulses. The black and blue (dark gray) 
solid lines correspond to the dynamics generated by the Lapert and the Krotov 
optimized pulses. The red (light gray) horizontal solid line depicts here the 
maximum permanent alignment that can be reached in the subspace T-ij, . (b) : 
X- Components of the optimized fields, (c): y- Components the optimized fields. 
The black and blue (dark gray) curves depict respectively the results derived 
with the Lapert and the Krotov algorithms. 



Acknowledgment 

The authors thank Daniel Reich for many helpful discussions and his comments 
on the manuscript. Financial supports from the Conseil Regional de Bourgogne 
and the QUAINT coordination action (EC FET-Open) are gratefully acknowl- 
edged. 

A Description of monotonically convergent algo- 
rithms 

The Krotov and Lapert optimization algorithms are summarized here for pure 
state quantum dynamics. This description is straightforwardly extended to 
the density matrix formalism. To simplify notation, we restrict ourself to the 
maximization of the projection onto a target state. The two algorithms can 
analogously be used for maximization or minimization of the expectation value 
of a given observable. The control problem is characterized by maximization of 
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the functional C, 

c = c, [{|v(t/))}] - r Ct [{|^(t))}, £;(i)] dt, (17) 

Jo 

where Ct , is the final time cost functional and Ct the running cost. The pa- 
rameter E denotes the external field and \il>{t)) the wave function describing the 
state of the system at time t. Its time evolution is governed by Eq. ^. If \il]f) 
is the target state, the final cost can be defined as follows: 

ct,[{i^(t/))}] = i(V'/iV'(i/))r (18) 

Here, only a running cost which does not depend on the state of the system is 
considered: 

Ctm{m,E{t)]=g[E{t)]. 

Extension to a state-dependent running cost is described in Ref |19) . The main 
difference between the two algorithms is in the choice of the running cost. 

A.l Krotov's method 

The derivation of the Krotov algorithm presented here follows closely Ref [TF] , 
specializing it to a non-linear interaction of the system with the control field, 
see Eqs. © and ([SJ. There is no requirement for a specific power of E{t) the 
running cost g so we choose it to minimize the change in the energy of the 
field [2T], 

9[E{t)]^^^{E{t)^Er.ff, (19) 

with Eref denoting a reference field, S(t) a shape function and A a weight. 

Krotov's method is based on the construction of an auxiliary functional, 
£ [{IV')}, -E, $] with $ an arbitrary functional. It is chosen such that the maxi- 
mization of C is equivalent to maximization of C of Eq (|T7l) . $ is used to ensure 
a global minimum with respect to changes in the state. Then any change in 
the state will lead to an increase in the value of C, i.e., to monotonic conver- 
gence [ini [2D] ■ This is achieved by expanding $ to second order in the change 
of the state, |A?A(i)), 

^[m}A = {x{t)mt)) + {i^{t)\x{t)) 

+ i(A^(t)k(i)|A^(t)). (20) 

When the equations of motion are linear with respect to the state and the run- 
ning cost functional does not depend on the state, the second order contribution 
is not required |19J. Since, in this work, we consider a quantum control problem 
which fulfills these conditions, a first order construction of $ is sufficient. The 
construction of $ is described in detail in Ref. [19]. The auxiliary functional is 
defined by: 

L[m},EM = G[{|V'(^/))}]-<i>[{|^(0))},0] 
ft, 

R[{\i,{t))},E{t),t]dt, (21) 

/o 
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where the final time functional and running functional, G and R, are given by 

Giimf))}) = a, [{i^(t/))}] + <i>[{i^(t/))},t/] , (22) 



+i(V'(t)|i/(V(.^l<i>) (23) 

For a maximization problem, the following conditions have to be fulfilled: 

C[{\^k)},Ek,'P] < C[{\iJk+i)},Ek+u'^], (24) 

where k indicates the iterative step. Sufficient conditions for maximizing C 
translate into maximizing G and minimizing R at each time: 



C[{\iJk+i)},Ek,^] - £[{m},Ek+i,^] = 

Ai+ / A2{t)dt+ / A3{t)dt, 
Jo Jo 

where the A^, i = 1, • • • , 3 are given by: 

A,^G[{\ijk+i{tf))}]-G[{\Mtf))}], 



(25) 



(26) 



and 



A2(t) = R[{\^k+i{m,Ekit),t] 

-R[{\^k+iim,Ek+i{t),t] 



A3(t) = R[{\Mm,Ek{t),t] 

-R[{\Mm,Ek+iit),t]. 



(27) 



(28) 



Non-negativeness of the A^ {i — 1, • • • , 3) ensures monotonic convergence. For 
a quantum control problem where the cost functional is state-independent and 
the equations of motion are linear with respect to the state, positivity of Ai and 
A3 is automatically satisfied PUI \n\ . Non-negativeness of A2 can be obtained 
by a proper choice of A and the shape function S{t) [19J: 



A . 

where Aj^^e is the spectral radius of Af^ == d^H [E{t)] /dE"^. 

Evaluating the extremum condition for C yields the control equations, 

1. Equation of the control field: 



(29) 



dE 



it) 



(30) 



Bfc4 



2am 



Xk{t) 



dH[E{t)] 
dE 



i^k+iit) 
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2. Equation of motion for the adjoint state |x), with 'initial' condition: 



i§i\Xk{t)) = H[E,{t)]\xk{t)) (31) 

\Xk{t^T)) = -V(^,|C*,||^^^. (32) 

3. Equation of motion of the state \^jJ), with initial condition IV'mi) 

z^|^,+i(i)) = H[Ek+i{t)]\i^k+i{t)) (33) 

|V'fe+i(i = 0)) = 1^,™). (34) 

At each iteration, the update of the field Ek+i{t) obtained from Eq. ([50)) involves 
a backward and a forward propagation, Eq. ([?T|) and Eq. (|55|) . respectively. If 
the system interacts non-linearly with the control field, both left and right hand 
sides of Eq. (j30l) depend on Ek+i{t). For simplicity, we assume that the change 
of the control field between iterations k and fc + 1 is small enough such that 
dH{Eit))/dE\E,^, « dH{E{t))/dE\E,. 

A. 2 The Lapert approach 

A. 2.1 General description 

While in the Krotov method, this running cost minimizes the change in the 
energy of the field, in the Lapert algorithm, this choice is different. Basically, g 
is chosen so as Eq. ([5n|) admits a real solution at any time t. The cost is defined 
as follows: 

g[E{t)]^^^{E{t)-Ereff\ (35) 

For a Hamiltonian given by Eq. ^, Eq. ([50)1 leads to 

2n^^{E,+,it)-E,it)f-'- 

23 [(xfclM + 2aEk+i + 3l3El^^\^k+i)] = (36) 

The left hand side of Eq. ([55)1 can be viewed as a polynomial of Ek+i- Choosing 
the integer n such that \/ S{t)E'^^^ is a monomial of order higher than the right 
hand side of Eq. psp ensures that there exists a real solution to the equation at 
each time t. For a non-linearity of order 3, n = 2 is sufficient. The conditions 
for monotonic convergence are determined through variation of AC given by: 

AC — Cfe+i — Cfe 

= Ct,mk+i{ts)))]~Ct,[{\Mtf))}] 
l-tf 

g[Ek+i{t)]-g[Ek{t)] dt. (37) 

which needs to be positive [25] . 
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Figure 8: (Color online) Variation of one of the real roots of Eq. (|55)) as a 
function of a; = {xit)\^{t)) for different values of A, 10^ < A < 10^. 



A. 2. 2 Role of the parameter A 

We discuss in this section the way the parameter A affects the optimized solution. 
For this purpose, we have analyzed the behavior of one of the real roots of the 
polynomial Eq. (|36p . To simplify the description, the operators fl, a and /3 have 
been replaced by their maximum eigenvalues. For n = 2, Eq. (j36p can then be 
written as follows: 



{Ek+i{t) - Ek{t)f - 6\(3\, 



S{t) 

-4|a|max3[(Xfc|V'fc+l>]^fe+l 
-2|/i|maxQ[(Xfc|V'fc+l>] =0. 



,^[{Xk\^k+i)]Ef 



fc+1 



(38) 



Figure [S] illustrates the variation of one of the real roots of Eq. ([55]) as a function 
of a; = {xit)\'4'{t)) for different values of A. The range of A is taken from 10^ to 
10^. For large values of A, the variation of the root is very slow with respect 
to X while for a value smaller then 10^, the change of the roots can be very 
fast. This observation qualitatively explains the fast oscillations occurring in 
the Lapert optimized fields. 
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